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We present direct numerical simulations (DNSs) of turbulent channel flow with passive Lagrangian 
polymers. To understand the polymer behavior we investigate the behavior of infinitesimal line 
elements and calculate, for the first time, the PDF of finite-time Lyapunov exponents and from 
them the corresponding Cramer's function for the channel flow. We study the statistics of polymer 
elongation for both the Oldroyd-B model (for Weissenberg number Wi < 1 ) and the FENE model. 
We use the location of the minima of the Cramer's function to define the Weissenberg number 
precisely such that we observe coil-stretch transition at Wi « 1. For linear polymers (Oldroyd-B 
model), with Wi < 1, we find that the PDF of polymer extensions shows power- law behavior with 
an exponent consistents with earlier analytical results pQ. For Wi > 1 (FENE model) the polymer 
are significantly more stretched near the wall than at the centre of the flow. Furthermore near the 
wall the polymers show a strong tendency to orient along the stream- wise direction of the flow but 
near the centerline the statistics of orientation of the polymers is consistent with analogous results 
obtained recently in homogeneous and isotropic flows [2]. 

PACS numbers: 



I. INTRODUCTION 

Turbulent flows with polymer additives have been an 
active field of interest since the early experiments of 
Toms [3 which showed the phenomenon of drag reduc- 
tion on the addition of small amounts (few parts per mil- 
lion) of long-chained polymers to turbulent wall-bounded 
flows. If disturbed from its equilibrium configuration it 
relaxes back and for simplest polymer models this relax- 
ation is governed by a single time scale r po i y . If such a 
polymer is then placed in a straining flow where the strain 
can be characterised by inverse of a time scale Tfl U id the 
polymer can go from its coiled state to a stretched state if 
the ratio of the two time scales, the Weissenberg number 
Wi > 1 [4 . Thus in turbulent flows with strong strain 
the polymers can go through coil-stretch transition; the 
stretched polymers can then make significant contribu- 
tion to the Reynolds stresses and this can result in drag- 
reduction [5j[6]. Thus to understand drag-reduction we 
must understand not only the mechanism of coil-stretch 
transition and but also the statistical properties of the 
stretched polymers. 

On the analytical side, the mechanism of coil-stretch 
transition in synthetic turbulent flows for simple dumbell 
models of polymers have been elucidated in Refs. [UEHH]. 
These studies use a flow that is white- in-time and smooth 
on the smallest scales - the so called Batchelor-Kraichnan 
flow - and linear and nonlinear models for polymers with 
a single characteristic relaxation time r po i y . Further they 
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ignore the backreaction of the polymers to the fluid. For 
linear polymers below coil-stretch transition Ref. [1] pre- 
dicts a power-law tail for the probability distribution 
function (PDF) of polymer extension. The exponent of 
the power-law a can be obtained via the Cramer's func- 
tion of the PDF of the finite-time Lyapunov exponent of 
Lagrangian particles in the flow, v.i.z., 

a = S"08 + — -//) (1) 

r poly 

where j3 must be obtained by solving the differential 
equation, 

Sifi + — -ft)- PS\P + — -p)=0 (2) 

r poly r poly 

Here S is the Cramer's function and ft is the location of 
the minima of S. These relationships have been numeri- 
cally verified for homogeneous and isotropic two dimen- 
sional flows [9] but not yet for wall bounded flows. 

In this paper we confirm the existence of this relation- 
ship by Direct Numerical Simulation (DNS) of a channel 
flow with polymer additives. The polymers are modelled 
by Lagrangian passive dumbells. We also calculate for 
the first time, the PDF of finite-time Lyapunov expo- 
nents for Lagrangian particles in a channel flow and the 
Cramer's function. We further find that the range of scal- 
ing shown by the PDF of polymer extensions depends on 
the wall-normal coordinate but the scaling exponent a is 
independent of the wall-normal coordinate. 

We use the location of the minima of the Cramer's 
function, /i, as the inverse characteristic time scale of the 
fluid to define our Weissenberg number as 

Wi /iTpoly (3) 
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Our simulations show coil-stretch transition for Wi ~ 1. 
For Wi > 1 it is not possible to obtain a stationary PDF 
for linear polymers. In this regime we use the nonlinear 
FENE (Finitely Extendable Nonlinear Elastic) model. 
As we use a Lagrangian model for the polymer we are 
able to use relatively high Wi. 

In addition we also find that the PDF of polymer ex- 
tensions depends on the wall-normal coordinate, v.i.z, the 
polymer are more stretched near the wall than at the cen- 
tre of the flow. We further study the orientation of the 
polymers with respect to the channel geometry and the 
local velocity gradient tensor. Our results show that the 
orientation of the polymers is predominantly determined 
by the inhomogeneity of the flow, i.e., the wall- normal 
coordinate as opposed to the local strain tensor. How- 
ever, for polymers near the centre of the channel we find 
that the orientation is also influenced by the principal di- 
rections of the rate-of-strain tensor, as has been seen in 
DNS of polymers in homogeneous isotropic flows [2j [10] . 

The rest of the paper is organised as follows. In the 
next section, Section [TT] we describe the equations we 
solve and the details of the numerical algorithm we use. 
Our results follow next in Section lllll which is divided into 
three parts. In Section [TlI A| we present our results on the 
PDF of finite-time-Lyapunov-exponents (FTLEs). The 
results described in this section are therefore indepen- 
dent of the polymer equation. The statistics of polymer 
extensions for the two models considered are presented 



in Section IIIB and IIIC The polymer orientation is 
characterized by calculating the correlations between the 
polymer end-to-end vector a nd flu id vorticity and the 
rate of strain tensor (Section |lIID ). We show that the 
polymers are oriented along the direction of the mean 
flow near the wall whereas near the center of the channel 
polymer orientation statistics is similar to that of homo- 
geneous, isotropic turbulence [2]. Finally we present our 
conclusions in Section HVl 



II. EQUATIONS AND NUMERICAL METHODS 

The fluid is described by the Navier-Stokes equations, 
d t u + u • Vu = vV 2 u + Vp (4) 



with the incompressibility constraint, 



0. 



(5) 



Here u is the fluid velocity, v the kinematic viscosity and 
p the pressure. We use no-slip boundary condition at 
the walls and periodic boundary condition at all other 
boundaries. We have chosen our units such that the con- 
stant density p = 1. The x axis of our coordinate system 
is taken along the stream-wise direction, the y axis along 
the wall-normal direction and z axis along the span-wise 
direction. For brevity, we shall often use the common 
notation, U = u X: V = u y and W = u z . The x, y and z 
dimensions of our channel are L x x L y x L z = 4tt x 2 x 27r, 
with resolution Nx x Ny x Nz = 128 x 129 x 128. 



The turbulent Reynolds number Re r = U*L/v = 180 
is defined by the friction velocity U* = ^/ov and L = 
Ly/2, the half-channel width, where 



dU. 

wall 

dy 



(6) 



is the shear stress at the wall [11 . In the following we 
non-dimensionalise velocity and distance by U + = U/U^ 
and y + = y/y* respectively, using the friction length 
y* = v/U*. Time will also be measured in the unit of 
the large-eddy-turnover-time tl = (/7 center /L) _1 , where 
Reenter j g avera g e streamwise velocity at the center of the 
channel. The large-scale Reynolds number defined by 
Re = UoL/u = 4200 where Uo is the centreline stream- 
wise velocity for the laminar flow of same mass flux. 

We solve Eqs. @ and Q by using the SIMSON [12] 
code which uses a pseudo-spectral method in space 
(Chebychev- Fourier). For time integration a third-order 
Runge-Kutta method is used for advection and uniform 
pressure gradient forcing term. The viscous term is dis- 
cretized using a Crank-Nicolson method [13]. In Fig. [I] 
we plot as a function of y + at the stationary state 

of our simulations, where (•) denotes averaging over the 
coordinate directions x and z and also over time. Fur- 
ther details about the code validation can be found in 
Ref. [H]. 




FIG. 1: (Colour online) Normalised mean stream- wise velocity 
(U + ) versus the wall- normal coordinate y + plotted for one 
half of the channel. 



We use a Lagrangian model for the polymers where we 
solve one stochastic differential equation (SDE) for each 
polymer molecule. This model uses several approxima- 
tions which are as follows [T5j[16]: (A) The centre-of-mass 
of a polymer molecule follows the path of a Lagrangian 
particle. (B) Even when fully stretched the polymer 
molecule is very small compared to the smallest scales of 
turbulence. This approximation is well justified [6 . (C) 
A polymer molecule is modelled by two beads separated 
by a vector which represents the end-to-end distance of 
the polymer molecule. (D) The forces acting on the beads 
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are stokes drag, restoring force of an overdamped spring 
with time scale r po i y , and thermal noise. To be specific, 
we track 7V p = 2.16 x 10 5 Lagrangian passive tracers in 
the flow by solving 



dtr*(t\t ,ri) =vi(t\t ,r J ). 



(7) 



Where ri(t\to,r J ) is the position of the j-th Lagrangian 
particle which was at position at time to and 
vi(t\to,r J ) is its velocity with j = 1,...,7V P . The La- 
grangian velocity of a particle, which is generally at 
an off- grid point, is obtained by tri-linear interpolation 
from Eulerian velocity at the neighboring grid points. 
Equation ^ is integrated by a third order Runge-Kutta 
scheme. Each of these Lagrangian particles represent a 
polymer molecule. For j-th Lagrangian particle the vec- 
tor representing the end-to-end distance is denoted by R) 
and obeys the following dynamical equation: 



dtRUt) 




Bl. 



(8) 



Here a 3 aj3 = dpvi(t\to, r J Q ) a: f(R^) is the restoring force 
of the polymer, r po i y is the characteristic decay time of 
the polymer and B ] is a Gaussian random noise with 
(B a ) = and (B a (t)B{3(t')) = 5 af3 5(t - t'). The pref- 
actor of the random noise is chosen such that in ab- 
sence of external flow, i.e., o^o = 0, the polymer at- 
tains thermal equilibrium, (W a R^) = R^d^/^t. Here 
(•) denotes averaging over the noise B. For the lin- 
ear Oldroyd-B model f(R) = -R/r poly . For the FENE 
model f(R) = -R/r poly {l - (R/R^) 2 }. Eq. (§ is 
also solved by a third order Runge-Kutta scheme except 
for the noise which is integrated by an Euler-Maruyama 
method [IT] . 

To compare with the analytical theory of Ref. [1 
we also need to calculate the PDF of finite-time Lya- 
punov exponents of Lagrangian particles in this flow. 
For this we need to calculate the rate at which two 
infinitismally separated Lagrangian particles diverge as 
time progresses. For this purpose we also calculate the 
evolution of an infinitesimal vector in our turbulent flow, 
given by the equations, 



dtSxl 



7 a, (3 ^^(3' 



(9) 



Where fix* is a vector carried by the j-th Lagrangian 
particle. This is of course the same equation obeyed by 
a Lagrangian polymer, Eq. (|8|, if the restoring force of 
the polymer and the Brownian noise are omitted. 

The correspondence between our Lagrangian descrip- 
tion and the Eulerian description of polymeric fluids is 
that in the latter the dynamical variable for the poly- 
mers is the symmetric positive definite (SPD) tensor 
C a p = (R a Rp). A DNS of the Eulerian description has 
certain difficulties [T8H2T] . Firstly the numerical schemes 
used must preserve the SPD nature of C a p. Secondly for 
high Weissenberg numbers large gradients of can 



develop which can lead to numerical instability of the 
code. Stability can generally be restored by employ- 
ing either shock-capturing schemes p~9j [2TH23] or by in- 
troducing dissipation in the Eulerian description of the 
polymer [24, 25 . Lagrangian methods [26j[27] are gener- 
ally able to avoid such numerical pitfalls and can attain 
higher Weissenberg number. On the other hand it is 
quite straightforward to incorporate the back-reaction of 
the polymer into the flow in the Eulerian model but is 
tricky in the Lagrangian model [27, 28 . Note finally that 
more complicated Lagrangian models have also been em- 
ployed where a single polymer is represented by a chain 
of beads connected by springs [2j [10] . 



III. RESULTS 
A. Finite-time Lyapunov Exponents 

How two infinitesimally separated Lagrangian particles 
diverge in a turbulent flow has been a central topic in 
turbulence research for a long time. See, e.g., Ref [29] for 
a recent review. Below we reproduce the essential points 
needed to apply such ideas to stretching of polymers in 
turbulence. For incompressible flows and for long time 
T, much larger than the correlation time of cr ] a ^, \5x* (T) | 
grows exponentially with time. We integrate Eq. ([9| for 
each Lagrangian particle over a finite time interval T and 
compute the the finite-time Lyapunov exponent (FTLE) 



In 



\5at(t)\ 



\5xi(t-T)\\ 



(10) 



For large T, T oo the PDF of FTLEs is conjectured 
to have a Large deviation form [TJ [30] [31] 



P(H T ) ~ exp[-TS( MT ) 



(11) 



where S(iht) is called the Cramer's function or the en- 
tropy function. The Cramer's function has earlier been 
calculated from direct numerical simulations of two [9]- 
and three [32] dimensional homogenous isotropic tur- 
bulence, turbulence in the presence of homogeneous 
shear [33] . and for hydromagnetic convection [34] . 

Below we present, for the first time, the analysis of the 
PDF of FTLEs for channel flows. 

Since the channel flow is not homogeneous in the wall- 
normal direction the statistics can, in principle, depend 
on y + . Hence we label our particles by their wall- normal 
coordinate (y + ) at the final position, i.e., at time T. 
While integrating the equations for 5x\ Eq. ([9]), we store 
the evolution of Sx^ and use this to calculate /i J T for each 
of SxK To calculate the PDF of /ix we gather statistics 
in two different ways. First we calculate the PDF of /ix 
for all particles at a fixed y + . Furthermore we run our 
simulations over several T and after each time interval T 
the particles are redistributed uniformly across the chan- 
nel and their initial separation vector 5x^ (t = 0) oriented 
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randomly. By definition then we generate a P(/ix,?/ + ) 
which depends on y + . The PDFs for two different values 
of one close to the wall, and one near the centerline, 
are respectively plotted in the left and the right panel 
of Fig. Q for several time intervals T. The peak and 
mean of the PDFs are always positive showing that it 
is more probable for \5x*\ to increase exponentially as a 
function of time. For small T the PDFs near the centre 
and the PDF near the wall are very different from each 
other. Significantly larger elongation is found for those 
elements that are located closer to the wall. However, the 
two PDFs approach each other for large T. This can also 
be seen by plotting the mean value of the PDFs for three 
different y+es as a function of time, Fig. The peak 
value also shows a similar trend, see inset in Fig. (J3|. 
Hence an unique Cramer's function independent of y + 
can be defined for the channel flow for only very large 
time when the PDFs for different y + merge with one an- 
other. In a channel flow the stress tensor depends 
strongly on the wall- normal coordinate. Thus for short 
T we can expect that the PDF of /ix depend of y + . Con- 
versely, when T becomes much larger than the typical 
time it takes for a particle to travel from a position near 
the wall to a position near the centreline, we expect /ix 
to be independent of y + . Let us call this typical time the 
exit time T ex i t . Surprisingly, we observe from our data 
that we need to have T ex i t > 80tl for /ix to be indepen- 
dent of T. An estimate of the time it takes for a particle 
to travel from the wall to the centre of the channel can 
be given by the ratio of the half-width of the channel to 
the friction velocity, Taction = (L y /2)/U* ~ 15 in our 
simulations. In units of this time T exit > 5Tf riction which 
provides a better estimate than tl- 

From PDF of /ix for large T we calculate the Cramer's 
function using Eq. 11 We normalise P(/ix) such that its 
integral over the range of /ix is unity. For T > T exit the 
Cramer's function S(fi) calculated at different times T 
are found to be independent of T as it should be. This is 
shown by the collapse of the Cramer's function calculated 
at different times for a fixed y + = 62 in Fig. [ZJ Further- 
more, the Cramer's function thus found is independent 
of y + . The connection between the Cramer's function 
and the PDF of end-to-end polymer distance was shown 
in Ref [1 for linear polymers and in Ref [7 for nonlinear 
polymers. We discuss such relations in the next section 
where it will turn out to be useful to have an algebraic 
expression for the Cramer's function. To obtain such an 
expression we fit a polynomial function of the form 

S(A0 = a 2 (/i - ft) 2 + a 3 (/i - ft) 3 + a 4 (/i - ft) 4 (12) 

to our numerical data for S(fi) averaged over all values 
of y + and extract the coefficients <22,a3, and above. 
To estimate the errors in the coefficients we use the 
same fit to S(ja) obtained for individual y + and quote 
the range of obtained from such fits as the error in a^. 
The best fit is also plotted in Fig. Q. The coefficients 
corresponding to the best fit and their errors are given in 
the caption of Fig. Q. 



B. Statistics of polymer extensions: Oldroyd-B 
model 



Before we present detailed results on statistics of poly- 
mer extension let us precisely define the Weissenberg 
number, Wi. In simulations the Weissenberg number is 
defined as the ratio of the characteristic time-scale of the 
polymer, r po i y over a characteristic time scale of the fluid. 
Different definitions of the characteristic time scale has 
been used in literature to define the Weissenberg num- 
ber. Refs. [2j[l0l[20] use the Kolmogorov time scale t v to 
define Weissenberg number. We denote this Weissenberg 
number by Wi^ = Tpoiy/r^ where is the Kolmogorov 
time scale. In this paper we principally use the following 
definition for Weissenberg number 



Wi /iTpoly 



(13) 



where ft is the location of the minima of the Cramer's 
function S(ja). Our choice has two principal advantages. 
Firstly in channel flows the Kolmogorov scale depends 
on the wall- normal coordinate and hence is not unique. 
Secondly and more importantly a proper choice of Weis- 
senberg number gives the coil-stretch transition of the 
polymer at Wi « 1 which is exactly what we obtain. To 
compare with earlier simulations, which were all done in 
homogeneous flows, we also calculate Wi^ a11 and Wi^ enter 
where we use the Kolmogorov time scale at the wall and 
at the centre of the flow respectively. We typically ob- 
tain, Wi™ an « 30Wi and Wi^ enter « 5Wi. The differ- 
ent values of Wi that we use are given below, in bracket 
we mention the corresponding values of Wi^ enter for 
easy comparison with earlier simulations of homogeneous 
and isotropic turbulence. For the Oldroyd-B model, 
Wi(Wi^ enter ) = 0.1(0.5), 0.2(1),0.3(1.5), and 0.5(2.5) 
and Wi(Wi^ enter ) = 0.1(0.5), 0.3(1.5),0.5(2.5), 1.5(7.5), 
2.5(12.5),3.5(17.5),4.5(22.5),5.5(27.5),7(35), and 10(50) 
for the FENE model. We use R = 10~ 7 , 10~ 8 , and 
#max/#o = 100 and 1000 for the FENE model. 

Let us first present the results for the Oldyroyd-B 
model. Here we expect to see a power-law behavior for 
the PDF of polymer extensions, Q(R) ~ R~ a ~ x [I] for 
large R. In general, the calculation of PDFs from numeri- 
cal data is plagued by errors originating from the binning 
of the data to make histograms. Thus it is often a diffi- 
cult task to extract exponents such as a from such PDFs. 
A reliable estimate of such an exponent can be obtained 
by using the rank-order method [35] to calculate the cor- 
responding cumulative probability distribution function, 



Q C (R) 



R 



(14) 



If the PDF has a scaling range the cumulative PDF also 
shows scaling, i.e., Q C (R) ~ R~ a . For different values of 
Wi these the cumulative PDFs for y + = 74 are plotted 
in Fig. [5] The cleanest power-law is seen for Wi = 0.5. 
So we choose this Weissenberg number for further de- 
tailed investigation. First we show that the exponent of 
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FIG. 2: (Colour online) (left panel) PDF of /x T near the wall (y + « 6) for several values of T, v.i.z, T = l.(A), 3.(«), 5.(D), 
35.(0), and 100. (— ). (right panel) PDF of {It near the centerline (y + ~ 180) for several values of T, v.i.z, T = l.(A), 3.(«), 
5.(D), 35.(0), and 100. (— ). Plots at some other intermediate values of T are consistent with this plot, but are not show here 
for clarity. All times are measured in the unit of tl. 




FIG. 3: (Color online) (/jlt) for different values of T for three 
different positions, in the channel, near the wall(«), near cen- 
terline (o) and at y + — 84 (a). All times are measured in the 
units of tl. 



the power-law (Wi = 0.5) a = 0.81 ± 0.02 does not de- 
pend on the y + although the range over which scaling 
is obtained does, Fig. [6] The exponent a is obtained by 
fitting a power-law for five different values of y + . The 
mean is reported as the exponent above and the stan- 
dard deviation from the mean is reported as the error. 

This exponent a can be obtained from the Cramer's 
function S(fi) using the set of couple equations Eq. ([I]) 
and Q [1 which we rewrite below, 



a = S'( 



r poly 



■fi) 



(15) 



FIG. 4: (Color online) The collapse of the Cramer's func- 
tions S(/jl) versus \i at y + — 62 for different times T, v.i.z, 
T = 35.(+),45.(o),55.(A),70.(D), and 100.(*). All times are 
measured in the units of tl. The continious line is the poly- 
nomial fit as given in Eq. [l2| with p, = 0.105[0.088 0.13] 
and a 2 = 3.55[3.09 4.35], a 3 = -12.60[-27.48 - 4.29] and 
a 4 = 39.64[3.84 90.07] 



equation, 

S(/?H 



r poly 



■J£)-PS'{P + 



r poly 



■p) = (16) 



where j3 must be obtained by solving the differential 



Had the Cramer's function been well approximate by a 
parabola of the form = (/x — /i) 2 /A, Eq. ([!]) would 
simplify to a — (2/A)(l/r po i y — p). We have checked that 
this quadratic approximation does not give accurate re- 
sult for a in our case. Using the algebraic expression for 
S given in Eq. 12 , we numerically solve Eqs. ([T]) and 
This give a = 0.9 ± 0.29 which agrees with the results 
obtained from the cumulative PDF of polymers within 
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FIG. 5: (Color online) Log- log plot of the cumuliative PDF 
Q C (R) of the polymer extensions R as a function of R for dif- 
ferent values of Wi, v.i.z, Wi = 0.05(«), 0.1(o), 0.2(D), 0.3(A), 
and 0.5(0). 




lative PDF we remove all the polymers for which R is so 
large that R/Ro > ^cutoff where we choose i? C utoff = 100 
and 1000. The resultant cumulative PDFs are plotted in 
Fig. [7] where the original cumulative PDF is also plotted 
for comparison. It can be seen that the scaling behavior, 
although present, is valid over a much smaller range. In 
the same figure we have also plotted the cumulative PDF 
for the FENE model with R m ^/Ro = 1000. This also 
shows scaling with a reduced range. Thus we expect that 
in experiments reminiscence of this scaling law should be 
visible although it may be difficult to detect because of 
reduced range of scaling. 




FIG. 7: (Color online) The cumulative PDF Q C (R) of the 
polymer extensions R as a function of R for the Oldroyd- 
B model (o), Oldroyd-B model with all polymers with 
R/Rmax > ^cutoff removed with i? C utoff = 10 4 («) and 
i^utoff = 100(D), FENE model with R ma */Ro = 10 3 (A) and 

With iCaxARo = 10 2 (A). 



FIG. 6: (Color online) Log- log plot of the cumuliative PDF 
Q C (R) of the polymer extensions R as a function of R for 
different y + , v.i.z, y+ = 8(o), 74(a), and 180(D). We fit a 
straight line to the data between the two dashed vertical lines 
to calculate a. This fit is shown as the black line. 



error bars. We note here that the a we calculate using 
the Cramer's function has large margin of error because 
the a depends sensitively on the coefficients a& in Eq. 12 
To find these coefficients accurately we need to know the 
Cramer's function accurately for a large range of its ar- 
gument not just the location of its minima. Numerically 
this is a difficult task and would require collecting data 
over very long times. 

Finally let us comment on the possible experimental 
determination of the exponent a. In practice no polymers 
are linear and in most cases the ratio of R max (maximum 
possible extension of the polymer) to Rq (the equilibrium 
length) ranges between 100 to 1000. To see the effect of a 
maximum extension, we first select one of the cumulative 
PDFs plotted in Fig.pl say for y + = 74. From this cumu- 



C. Statistics of polymer extensions: FENE model 

Till now we have described the polymer statistics for 
Wi < 0.5. As we increase the Wi and make it close to 
unity no stationary statistics of the polymers is obtained. 
We interpret this by noting that we are close to the coil- 
stretch transition. A stationary state can be obtained 
by either including the feedback for the polymers into 
the fluid or by using nonlinear polymers e.g., the FENE 
model. We choose the second option. In the FENE model 
we have used R max /Ro = 100 and 1000. Our results as 
reported below does not depend on this parameter. 

First we show the coil-stretch transition by plotting 
the PDF of polymer extensions for three different Weis- 
senberg numbers Wi = 0.5, 1.5 and 10 in Fig. (8|. The 
sequence of these figures illustrate the coil-stretch transi- 
tion clearly. In the first plot the PDF is strongly peaked 
near zero which corresponds to the coiled state. The asso- 
ciated cumulative PDF is displayed in Fig. [7] and shows 
power-law behavior over a limited range of extensions. 
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The data for Wi = 1.5 are shown in the second figure: 
although the peak of the PDF is still close to zero the 
PDF is well spread over the whole range. The PDF for 
polymers with Wi = 10 has peak near R m &x which is 
the stretched state of the polymer. In these three fig- 
ures the PDF is plotted for y+ = 74 but the PDF at 
other wall-normal coordinates in the channel shows the 
same qualitative nature as was seen in the case of the 
Oldroyd-B model in Fig. 

In Fig. [9] we show how the mean polymer extension 
(R)xzi where the averaging is over the stream- wise and 
the span- wise direction, changes with Wi across the chan- 
nel for R max /Ro = 100. For a given Wi the average 
polymer extension is small near the wall, increases to a 
maximum around y + ~ 10 (this corresponds to the re- 
gion of maximum strain), and then decreases towards 
the center of the domain where the flow is close to ho- 
mogeneous turbulence. A similar trend is also seen for 
^max/^o = 1000. This trend has been seen in earlier 
DNS of polymeric turbulence in channel flows, see, e.g., 
Ref. [36] and the references therein. But for larger values 
of Wi the average polymer extension becomes almost uni- 
form across the channel (except very near the wall where 
it is always small). This is because the polymers that are 
stretched close to the wall on reaching the centerline are 
not able to relax fast enough as the polymer relaxation 
time scales are much larger than the fluid time scales. 




max 



FIG. 8: (Color online) The PDF of polymer extensions Q(R) 
as a function of R for different Wi showing the coil-stretch 
transition. The line with (x) symbols is for Wi = 0.5 (T po i y = 
5), the continuous line is for Wi = 10 (T po i y = 100) and the 
inset is for Wi = 1.5 (r po i y = 15). The PDF for Wi = 10 is 
multiplied by 2 to make it clear in the same scale. The dashed 
line shows power-law scaling with exponent a = 1.48. 
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FIG. 9: (Color online) The average polymer exten- 
sions (R)xz as a function of the wall-normal coordi- 
nate y + for different Weissenberg numbers, v.i.z., Wi = 
1.5(*), 2.5(B), 3.5(D), 4.5(A), 5.5(a), 7(«), 10(o). The maxi- 



mum occurs at 



10. 




FIG. 10: (Color online) The maxima of the average polymer 
extensions {R) X z across the channel as a function of Wi in 
semilog scale. 



D. Statistics of polymer orientation 

In this section we present the results related to the 
orientation of the polymers. First let us discuss the ori- 
entation of the polymers with respect to the geometry 
of the channel. Let us denote the unit vector along R 
to be e. The PDF of the three components of e, e x , 
e y and e z , (i.e., three direction cosines of R) are plot- 
ted in the left panel of Fig. (11) for polymers close to 
the wall (y + « 7) and for three different values of Wi. 
For Wi < 1 , i.e., below the coil-stretch transition the 



polymers are almost equally probable to point in any di- 
rection or in other words as the polymers are coiled as a 
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sphere no preferential direction is picked up. Above the 
coil-stretch transition polymers close to the wall have a 
high probability of being oriented along the x axis, which 
is the stream-wise direction. A similar plot for polymers 
close to the centreline (y + ~ 180) is given in the right 
panel of Fig. ( 11 ). For small Wi all directions are equally 



probable. But as Wi increases here too the polymers get 
preferentially oriented along the stream-wise direction al- 
though the trend is much weaker than what is seen near 
the wall. 



We have also investigated the orientation of the poly- 
mers with respect to the three principal directions of the 
rate of strain tensor. For this purpose we first determine 
the three real eigenvalues of the symmetric rate of strain 
tensor and order them such that Ai > A2 > A3. We de- 
note the components of the unit vector e (which is the 
unit vector along R) along these three perpendicular di- 
rections by ei, e2 and e3, these are merely the cosines of 
the angles between R and the three principal directions 
of the strain tensor. The PDFs of ei, 62 and es are plot- 
ted in the left panel of Fig. (12) for polymers close to the 
wall (y + ~ 7) and for three different values of Wi. The 



peak seen in Fig. (12) corresponds to the polymers orien- 



tating along the stream-wise direction as shown already 
in Fig. (11). Interestingly the polymers are not preferen- 



tially oriented along the strongest direction of strain Ai 
but along the stream- wise direction. This has an angle 
of about 45 degrees with respect to the x-axis since the 
main component of the strain rate is because of the wall- 
normal shear dU/dy. Close to the centerline however 
the PDFs look quite different [right panel of Fig. (12)]. 
For small Wi there is no preferential orientation but as 
Wi increases the polymers develops a trend of orienting 
parallel to the direction of either Ai or A2 and shows 
ant i- alignment to A3. 



Finally we look at the relative orientation between the 
polymer end-to-end vector R and the vorticity vector uj. 
Close to the wall we find that PDF of the cosine of the 
angle ijj between R and uj has a peak at zero, see left 
panel of Fig. (13). This implies that the polymers show 



a weak tendency to lie in the plane perpendicular to u;. 
However this trend is reversed near the centerline [right 



panel of Fig. (13)] where the polymers orient along the 
vorticity vector. 



IV. CONCLUSIONS 



We have presented in this paper an extensive study 
of the passive Lagrangian polymers in turbulent channel 
flow by Direct Numerical Simulations. We have used 
both linear (Oldroyd-B) and nonlinear (FENE) poly- 
mers. To understand the statistics of polymer end-to-end 
vector it is necessary to know the statistics of the Finite 
Time Lyapunov Exponents. For this purpose in addition 
to the polymers we have solved the equation of evolution 
of infinitismal line elements in the turbulent flow and 
for the first time calculated the FTLEs for an inhomo- 
geneous flow. We find that at very large time the PDF 
of FTLEs admit a large deviation description and a cor- 
responding Cramer's function can be found. In addition 
we use the location of the minima of the Cramer's func- 
tion to define our Weissenberg number. Consequently 
for the FENE model we observe coil-stretch transition 
at Wi ~ 1. For the Oldroyd-B model we find that the 
PDF of polymer extension shows power-law behavior for 
Wi < 1. We calculate the exponent of this power-law us- 
ing the rank-order method. We also calculate the same 
exponent from the Cramer's function using the theory of 
Ref. [l] . These two different calculations match within er- 
ror, validating the theory of Ref pQ. This shows that the 
idealizations used in Ref. pQ , in particular the assumption 
that in Lagrangian coordinates the rate-of-strain tensor 
a a i3 is delta-correlated in time is a reasonable approxi- 
mation at least for linear polymers below the coil-stretch 
transition even in the case of a realistic flow. For the 
FENE model we cannot meaningfully calculate the PDF 
of polymer extension from the Cramer's function using 
the results of Ref. [7 because our numerically calculated 
Cramer's function is not accurate enough for this exer- 
cise. For the FENE model we find that the polymers 
are more extended near the wall, but the difference de- 
creases as Weissenberg number increases far beyond the 
coil-stretch transition. We further find that near the cen- 
tre of the channel the orientational statistics of the poly- 
mers show similarity to orientational statistics obtained 
for homogeneous and isotropic flows [2], i.e., they align 
along either of the two largest directions of strain and 
tend to orient orthogonal to the third principal direction 
of strain. But a much stronger orientational trend is seen 
near the wall where the orientations of the polymers are 
along by the stream-wise direction. 



To summarise the polymers near the wall shows 
the cleanest trend in their orientation. They show a 
strong tendency to line along the stream- wise directions. 
Weaker trends are seen near the center. The statistics 
of orientation of polymers near the center of our flow is 
very similar to the statistics of orientation of polymers 
obtained in homogeneous and isotropic flows [2 j. But 
note that the orientation effects we see near the wall are 
much stronger compared to those near the centerline. 



Although our DNSs involve passive polymers it is 
possible to have insights on polymeric drag reductions 
from these simulations. We can calculate the polymeric 
stress from our simulations and add this to the Reynolds 
stresses to see how they change the Reynolds equations 
of the flow. It would be interesting to see how much of 
drag-reduction can be described by this simple approach. 
Such results will be presented in a future publication. 
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FIG. 11: (Color online) PDF of the three direction cosines of polymer end-to-end separation vector R for polymers near the 
wall (left panel) and for polymers at the center of the channel (right channel). Three different values of Wi are used. Namely, 
Wi = 0.1(o), 1.5(A), 4.5(B). The data for Wi = 1.5 and 4.5 coincide on each other. The PDFs of e x and e y are respectively 
plotted using continuous line with symbols (P x ) and dashed lines with symbols (P y ). The inset shows the PDF of e z , P z . 




FIG. 12: (Color online) PDF of ei, e2 and e3, components of the unit vector along R along the three principal directions 
of strain, for polymers near the wall (left panel), and for polymers near the centerline (right panel). Three different Wi are 
used. Namely, Wi = 0.1(o), 1.5(A), 4.5(B). The data for Wi = 1.5 and 4.5 coincide on each other. The PDFs of e± and e<2 are 
respectively plotted using continuous line with symbols (Pi) and dashed line with symbols (P2). The inset shows the PDF of 
e 3 , Ps. 
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